#!/usr/local/bin/python2.6

from ruffus import *
import os
import sys
import re

exe_path = os.path.split(os.path.abspath(sys.argv[0]))[0]
sys.path.insert(0, os.path.abspath(os.path.join(exe_path,"..", "..","..")))

GATK="java -Xmx11g -jar /scratch/FA_test/GATK/GenomeAnalysisTK.jar "
REF="/scratch/FA_test/data/human37.fasta"
DBSNP="--DBSNP /scratch/FA_test/data/snp131.dat"
DBSNP2="-B:dbsnp,AnnotatorInputTable /scratch/FA_test/data/snp131.dat"
AL1000="-J 1000g,/scratch/FA_test/data/1000genomes_allelic_freq.dat,1000g.ID=dbsnp.name"
REFSEQ="-B:refseq,AnnotatorInputTable /scratch/FA_test/data/refGene-big-table-hg19.txt"
REFL="-J reflink,/scratch/FA_test/data/refLink.hg19.txt,reflink.mrnaAcc=refseq.name"
CCDS="-B:ccds,AnnotatorInputTable /scratch/FA_test/data/ccdsGene_protein.txt"
OUTF="-s dbsnp.refUCSC,dbsnp.observed,dbsnp.strand,dbsnp.name,dbsnp.func,dbsnp.valid,1000g.allelicFr,phC.name,phC.score,reflink.omimId,refseq.name,refseq.name2,refseq.positionType"
PH46="-B:phC,AnnotatorInputTable /scratch/FA_test/data/phastCons.dat"



@transform("*.bam", suffix(".bam"), ".recall.csv")
def gatk_covariance(bamFile, recFile):
    os.system("%s -R %s -I %s -T CountCovariates %s -cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate  -recalFile %s" %(GATK, REF, bamFile, DBSNP, recFile))
        
@follows(gatk_covariance)
@transform("*.bam", suffix(".bam"), ".recal.bam")
def gatk_recalibrate(bamFile, outFile):
    bf = re.split('\.', bamFile)
    recFile = bf[0] + ".recall.csv"
    os.system("%s -R %s -I %s -T TableRecalibration --out %s -recalFile %s" % (GATK, REF, bamFile, outFile, recFile))
    os.remove(recFile)

@follows(gatk_recalibrate)
@transform("*.recal.bam", suffix(".recal.bam"), ".realign.bam")
def gatk_realign(bamFile, realignBam):
    tmpFile = "bam_recal.list"
    bf = re.split('\.', bamFile)
    sortBam = bf[0] + ".recal.sorted.bam"
    tmpSort = bf[0] + ".recal.sorted"
    os.system("samtools sort %s %s" % (bamFile, tmpSort))
    os.system("samtools index %s" % (sortBam))
    os.system("%s -T RealignerTargetCreator -I %s -R %s -o %s" % (GATK, sortBam, REF, tmpFile))
    os.system("%s -I %s -R %s -T IndelRealigner -targetIntervals %s --out %s" % (GATK, sortBam, REF, tmpFile, realignBam))
#    os.remove(bamFile)
#    os.remove(tmpSort)
#    os.remove(tmpFile)
#    os.remove(sortBam)


@follows(gatk_realign)
@transform("*.realign.bam", suffix(".realign.bam"), ".vcf")
def gatk_snp_call(bamFile, vcfFile):
    bf = re.split('\.', bamFile)
    sortBam = bf[0] + ".sorted.bam"
    tmpSort = bf[0] + ".sorted"
    os.system("samtools sort %s %s" % (bamFile, tmpSort))
    os.system("samtools index %s" % (sortBam))            
    os.system("%s -R %s -T UnifiedGenotyper -baq CALCULATE_AS_NECESSARY -I %s --out %s -stand_call_conf 30.0" % (GATK, REF, sortBam, vcfFile))
    os.remove(bamFile)
    
@follows(gatk_snp_call)
@transform("*.vcf", suffix(".vcf"), ".vcf.annot")
def gatk_annotate(vcfFile, annotFile):
    os.system("%s -T GenomicAnnotator -R %s -B:variant,vcf %s %s %s %s %s %s %s -o %s -BTI variant -maxJoin 16000000 %s" % (GATK, REF, vcfFile, DBSNP2, AL1000, REFSEQ, REFL, CCDS, PH46, annotFile, OUTF))

pipeline_run([gatk_snp_call])
